#### Replication Code for:
#### Ministries Matter: Technocrats and regime loyalty under autocracy
#### Author: Erin York
#### Created: 2/3/2023


rm(list = ls())

# required packages
library(tidyverse)
library(MatchIt)
library(sandwich)
library(survival)
library(stargazer)
library(lubridate)
library(broom)
library(logr)

options("logr.autolog" = TRUE)

# Open the log
lf <- log_open(file.path("logmain.log"))

# Send message to log
log_code()


# load in data ------------------------------------------------------------

# set working directory to location of query_data.csv
df<- read.csv("query_data.csv")
df<- df %>%
  mutate(exactpost = ymd(exactpost))

# date of cabinet switch:
minchange<- ymd("2013-10-10")


# exact matched analysis - MAIN -------------------------------------------

matchdat<- df %>%
  # remove data missing party data
  filter(!is.na(party_group)) %>%
  # remove non-partisan/nontechnocrat ministers
  filter(!(party_min ==0 & techno_min ==0))

## exact matching for techno/partisan queries:
dfregex_main<- matchit(techno_min ~ # timing attributes
                         as.factor(year) +late_in_term + 
                         # ministry/cabinet attributes
                        first_cab+ sovereignty + large_min +
                         # query content attributes
                         georeference + casework  +critical + 
                         # asker attributes
                         asker_leader + asker_natlist + factor(party_group), 
                       method = "exact", 
                       estimand = "ATT",
                       data = matchdat)

dta_regex_main <- match.data(dfregex_main) %>%
  left_join(df)

dim(dta_regex_main)
summary(dfregex_main)
# 93% of treated
# 71% of untreated


## Table 1 specifications - OLS
t1_m1<- (lm(received_resp ~ techno_min , dta_regex_main, weights = weights))

t1_m2<- (lm(received_resp ~ techno_min + 
              # timing attributes
              as.factor(year) +late_in_term + 
              # ministry/cabinet attributes
             first_cab+ sovereignty + large_min +
              # query content attributes
              georeference + casework  +critical + 
              # asker attributes
              asker_leader + asker_natlist + party_group, 
            dta_regex_main, weights = weights))

# year FE + other attributes
t1_m3<- (lm(received_resp ~ techno_min +  mincode +
              as.factor(year) +late_in_term + # timing
             first_cab+ # ministry/cabinet
              georeference + casework  +critical + # content
              asker_leader + party_group +asker_natlist, # asker
            dta_regex_main, weights = weights))

# ses clustered at min-period level
se_t1_m1<- sqrt(diag(vcovCL(t1_m1, cluster = ~ min_period)))
se_t1_m2<- sqrt(diag(vcovCL(t1_m2, cluster = ~ min_period)))
se_t1_m3<- sqrt(diag(vcovCL(t1_m3, cluster = ~ min_period)))



# survival analysis

dta_regex_cox <- dta_regex_main %>% 
  # filter out 48 observations with incorrect response date (life_pre <= 0)
  filter(life_pre >= 0)

## estimate survival models
mod<- Surv(dta_regex_cox$life_pre, dta_regex_cox$received_resp)

t1_m4<- (coxph(mod ~ techno_min, dta_regex_cox, cluster = mincode, weights = weights))

t1_m5<- coxph(mod ~ techno_min + as.factor(year) +late_in_term + # timing
               first_cab+ sovereignty + large_min +# ministry/cabinet
                georeference + casework  +critical + # content
                asker_leader + party_group +asker_natlist  , # asker
              dta_regex_cox, 
              cluster = min_period, weights = weights)


t1_m6<- coxph(mod ~ techno_min + mincode + as.factor(year) +late_in_term + # timing
               first_cab+ # ministry/cabinet
                georeference + casework  +critical + # content
                party_group + asker_leader  + asker_natlist, # asker
              dta_regex_cox, cluster = min_period, weights = weights)

## TABLE 1 ##
sink("Table1.tex")
stargazer(t1_m1, t1_m2, t1_m3, t1_m4, t1_m5, t1_m6,
          se = list(se_t1_m1, se_t1_m2, se_t1_m3,
                    summary(t1_m4)$coefficients[,4], 
                    summary(t1_m5)$coefficients[,4],
                    summary(t1_m6)$coefficients[,4]),
          keep = "techno*",
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("IPW:", rep("\\checkmark", 6)),
                           c("Controls:", rep( c("", "\\checkmark", "\\checkmark"),2)),
                           c("Ministry FE:", rep( c("", "", "\\checkmark"),2)),
                           c("Model:", rep("OLS", 3), rep("Cox PH", 3))),
          model.names = F)
sink()
####


## TABLE A4 ##
sink("TableA4.tex")
stargazer(t1_m1, t1_m2, t1_m3, t1_m4, t1_m5, t1_m6,
          se = list(se_t1_m1, se_t1_m2, se_t1_m3,
                    summary(t1_m4)$coefficients[,4], 
                    summary(t1_m5)$coefficients[,4],
                    summary(t1_m6)$coefficients[,4]),
          omit = c("mincode*", "party*", "year*", "x"),
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister",
                               "Late in term", 
                               "2012 cabinet",
                               "Ministry of sovereignty",
                               "Large ministry",
                               "Georeference", "Casework", "Critical",
                               "Asker - leader", "Asker - national list"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c( "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("IPW:", rep("\\checkmark", 6)),
                           c("Year FE:", rep( c("", "\\checkmark", "\\checkmark"),2)),
                           c("Party FE:", rep( c("", "\\checkmark", "\\checkmark"),2)),
                           c("Ministry FE:", rep( c("", "", "\\checkmark"),2)),
                           c("Model:", rep("OLS", 3), rep("Cox PH", 3))),
          model.names = F,
          no.space = T)
sink()
####




# matched mechanism analysis - MAIN ---------------------------------------

### subset based on georeferences

## models including only queries with a georeference
geo_m1<- (lm(received_resp ~ techno_min , 
             dta_regex_main %>% filter(georeference == 1|casework == 1), weights = weights))

geo_m2<- (lm(received_resp ~ techno_min + as.factor(year) +late_in_term + # timing
              first_cab+ sovereignty + large_min +# ministry/cabinet
               critical + # content
               party_group + asker_leader  + asker_natlist , # asker
             dta_regex_main %>% filter(georeference == 1|casework == 1), weights = weights))

# ses clustered at min-period level
se_geo_m1<- sqrt(diag(vcovCL(geo_m1, cluster = ~ min_period)))
se_geo_m2<- sqrt(diag(vcovCL(geo_m2, cluster = ~ min_period)))

# models including only queries without georeferences
geo_m3<- (lm(received_resp ~ techno_min , 
             dta_regex_main %>% filter(georeference == 0&casework == 0), weights = weights))

geo_m4<- (lm(received_resp ~ techno_min + as.factor(year) +late_in_term + # timing
              first_cab+ sovereignty + large_min +# ministry/cabinet
               critical + # content
               party_group + asker_leader  + asker_natlist , # asker
             dta_regex_main %>% filter(georeference == 0&casework == 0), weights = weights))

# ses clustered at min-period level
se_geo_m3<- sqrt(diag(vcovCL(geo_m3, cluster = ~ min_period)))
se_geo_m4<- sqrt(diag(vcovCL(geo_m4, cluster = ~ min_period)))

## TABLE 3 ##
sink("Table3.tex")
stargazer(geo_m1, geo_m2, 
          geo_m3, geo_m4,
          se = list(se_geo_m1, se_geo_m2,
                    se_geo_m3, se_geo_m4),
          keep = "techno*",
          dep.var.labels = c("Received Response"),
          covariate.labels = c("Technocrat Minister"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("IPW:", rep("\\checkmark", 4)),
                           c("Controls:", rep( c("", "\\checkmark"),2)),
                           c("Subset:", rep("Const.", 2), rep("Non-const.", 2))
          ),
          model.names = F)
sink()
####

## subset based on sensitive/critical queries

## models including only sensitive/critical queries
crit_m1<- (lm(received_resp ~ techno_min , 
              dta_regex_main %>% filter(critical == 1), weights = weights))

crit_m2<- (lm(received_resp ~ techno_min + as.factor(year) +late_in_term + # timing
               first_cab+ sovereignty + large_min +# ministry/cabinet
                georeference + casework  +critical + # content
                party_group + asker_leader  + asker_natlist , # asker
              dta_regex_main %>% filter(critical == 1), weights = weights))

# ses clustered at min-period level
se_crit_m1<- sqrt(diag(vcovCL(crit_m1, cluster = ~ min_period)))
se_crit_m2<- sqrt(diag(vcovCL(crit_m2, cluster = ~ min_period)))


# ## models including only queries without sensitive/critical language
crit_m3<- (lm(received_resp ~ techno_min , 
              dta_regex_main %>% filter(critical == 0), weights = weights))

crit_m4<- (lm(received_resp ~ techno_min + as.factor(year) +late_in_term + # timing
               first_cab+ sovereignty + large_min +# ministry/cabinet
                georeference + casework  +critical + # content
                party_group + asker_leader  + asker_natlist , # asker
              dta_regex_main %>% filter(critical == 0), weights = weights))

# ses clustered at min-period level
se_crit_m3<- sqrt(diag(vcovCL(crit_m3, cluster = ~ min_period)))
se_crit_m4<- sqrt(diag(vcovCL(crit_m4, cluster = ~ min_period)))


## TABLE 4 ##
sink("Table4.tex")
stargazer(crit_m1, crit_m2, 
          crit_m3, crit_m4, 
          se = list(se_crit_m1, se_crit_m2,
                    se_crit_m3, se_crit_m4),
          keep = "techno*",
          dep.var.labels = c("Received Response"),
          covariate.labels = c("Technocrat Minister"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("IPW:", rep("\\checkmark", 6)),
                           c("Controls:", rep( c("", "\\checkmark"),2)),
                           c("Subset:", rep("Critical", 2), rep("Not Critical", 2))),
          model.names = F)
sink()
####



# matched analysis - ROBUSTNESS -------------------------------------------

### first robustness check - nearest neighbor match ##
set.seed(123456)
dfregex_rob1<- matchit(techno_min ~ # timing attributes
                         as.factor(year) +late_in_term + 
                         # ministry/cabinet attributes
                        first_cab+ sovereignty + large_min +
                         # query content attributes
                         georeference + casework  +critical + 
                         # asker attributes
                         asker_leader + asker_natlist + factor(party_group), 
                       method = "nearest", data = matchdat)


dta_regex_rob1 <- match.data(dfregex_rob1) %>%
  left_join(df)
summary(dfregex_rob1)

# OLS specifications with and without controls
rob_m1<- (lm(received_resp ~ techno_min , dta_regex_rob1))

rob_m2<- (lm(received_resp ~ techno_min + 
               as.factor(year) +late_in_term + # timing
              first_cab+ sovereignty + large_min +# ministry/cabinet
               georeference + casework  +critical + # content
               asker_leader + party_group +asker_natlist  , # asker
             dta_regex_rob1))

se_rob_m1<- sqrt(diag(vcovCL(rob_m1, cluster = ~ min_period)))
se_rob_m2<- sqrt(diag(vcovCL(rob_m2, cluster = ~ min_period)))


### second robustness check - alternative "sovereignty" measure (sovereignty_alt) ##
dfregex_rob2<- matchit(techno_min ~ # timing attributes
                         as.factor(year) +late_in_term + 
                         # ministry/cabinet attributes
                        first_cab+ sovereignty_alt + large_min +
                         # query content attributes
                         georeference + casework  +critical + 
                         # asker attributes
                         asker_leader + asker_natlist + factor(party_group),
                       method = "exact", data = matchdat)


dta_regex_rob2 <- match.data(dfregex_rob2) %>%
  left_join(df)

dim(dta_regex_main)
summary(dfregex_rob2)
# 96% of treated
# 23% of untreated (no ministries that regime didn't hold)

# OLS specifications with and without controls
rob_m3<- (lm(received_resp ~ techno_min , dta_regex_rob2, weights = weights))

rob_m4<- (lm(received_resp ~ techno_min + as.factor(year) +late_in_term + # timing
              first_cab+ sovereignty_alt + large_min +# ministry/cabinet
               georeference + casework  +critical + # content
               party_group + asker_leader  + asker_natlist , # asker
             dta_regex_rob2, weights = weights))

se_rob_m3<- sqrt(diag(vcovCL(rob_m3, cluster = ~ min_period)))
se_rob_m4<- sqrt(diag(vcovCL(rob_m4, cluster = ~ min_period)))



### third robustness check - fewer controls/larger dataset ##
dfregex_rob4<- matchit(techno_min ~ 
                         # timing attributes
                         late_in_term + 
                         # ministry/cabinet attributes
                        first_cab+ sovereignty + large_min +
                         # content attributes
                         georeference + casework  +critical +
                         # asker attributes
                         asker_leader  + asker_natlist , 
                       method = "exact", data = matchdat)


dta_regex_rob4 <- match.data(dfregex_rob4) %>%
  left_join(df)

dim(dta_regex_main)
summary(dfregex_rob4)
# 99% of treated
# 79% of untreated

# OLS specifications with and without controls
rob_m5<- (lm(received_resp ~ techno_min , dta_regex_rob4, weights = weights))

rob_m6<- (lm(received_resp ~ techno_min +
               late_in_term + # timing
              first_cab+ sovereignty + large_min +# ministry/cabinet
               georeference + casework  +critical + # content
               asker_leader  + asker_natlist , dta_regex_rob4, weights = weights))

se_rob_m5<- sqrt(diag(vcovCL(rob_m5, cluster = ~min_period)))
se_rob_m6<- sqrt(diag(vcovCL(rob_m6, cluster = ~min_period)))

## TABLE A6 ##
sink("TableA6.tex")
stargazer(rob_m1, rob_m2, 
          rob_m3, rob_m4, 
          rob_m5, rob_m6,
          se = list(se_rob_m1, se_rob_m2, 
                    se_rob_m3, se_rob_m4,
                    se_rob_m5, se_rob_m6),
          keep = "techno*",
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("Matching:", rep("Nearest", 2), rep("Exact Alt 1", 2), rep("Exact Alt 2", 2)),
                           c("IPW:", rep("", 2), rep("\\checkmark", 4)),
                           
                           c("Controls:", rep( c("", "\\checkmark"),3)),
                           c("Model:", rep("OLS", 6))),
          model.names = F)
sink()
####



## robustness - excluding co-partisan queries
## drop party identifier 
dfregex_noalign<- matchit(techno_min ~ as.factor(year) +late_in_term + # timing
                           first_cab+ sovereignty + large_min +# ministry/cabinet
                            georeference + casework  +critical + # content
                            asker_leader  +asker_natlist +
                            copartisan, # asker
                          method = "exact", data = matchdat)

dta_regex_noalign <- match.data(dfregex_noalign) %>%
  left_join(df)

summary(dfregex_noalign)

# OLS models
noal_m1<- (lm(received_resp ~ techno_min , dta_regex_noalign, weights = weights))


noal_m2<- (lm(received_resp ~ techno_min +
                as.factor(year) +late_in_term + # timing
               first_cab+  sovereignty + large_min +# ministry/cabinet
                georeference + casework  +critical + # content
                asker_leader  + asker_natlist, # asker
              dta_regex_noalign, weights = weights))

noal_m3<- (lm(received_resp ~ techno_min +mincode +
                as.factor(year) +late_in_term + # timing
               first_cab+ # ministry/cabinet
                georeference + casework  +critical + # content
                party_group + asker_leader  + asker_natlist, # asker
              dta_regex_noalign, weights = weights))

# ses clustered at min-period level
senoalign_m1<- sqrt(diag(vcovCL(noal_m1, cluster = ~ min_period)))
senoalign_m2<- sqrt(diag(vcovCL(noal_m2, cluster = ~ min_period)))
senoalign_m3<- sqrt(diag(vcovCL(noal_m3, cluster = ~min_period )))


# survival

dta_regex_cox2 <- dta_regex_noalign %>% mutate(x = 1) %>% filter(life_pre >= 0)
mod<- Surv(dta_regex_cox2$life_pre, dta_regex_cox2$received_resp)

noal_m4<- (coxph(mod ~ techno_min + x, dta_regex_cox2, 
                 cluster = min_period, weights = weights))


noal_m5<- coxph(mod ~ techno_min +  as.factor(year) +late_in_term + # timing
                 first_cab+ sovereignty + large_min +# ministry/cabinet
                  georeference + casework  +critical + # content
                  asker_leader  + asker_natlist, # asker
                dta_regex_cox2,  cluster = min_period, weights = weights)


noal_m6<- coxph(mod ~ techno_min +mincode+ as.factor(year) +late_in_term + # timing
                 first_cab+ sovereignty + large_min +# ministry/cabinet
                  georeference + casework  +critical + # content
                  asker_leader  + asker_natlist, # asker
                dta_regex_cox2,  cluster = min_period, weights = weights)



## TABLE A9 ##
sink("TableA9.tex")
stargazer(noal_m1, noal_m2, noal_m3, 
          noal_m4, noal_m5, noal_m6,
          se = list(senoalign_m1, senoalign_m2, senoalign_m3,
                    summary(noal_m4)$coefficients[,4],
                    summary(noal_m5)$coefficients[,4],
                    summary(noal_m6)$coefficients[,4]),
          keep = "techno*",
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("IPW:", rep("\\checkmark", 6)),
                           c("Controls:", rep( c("", "\\checkmark", "\\checkmark"),2)),
                           c("Ministry FE:", rep( c("", "", "\\checkmark"),2)),
                           c("Model:", rep("OLS", 3), rep("Cox PH", 3))),
          model.names = F)
sink()
####


# DiD analysis - MAIN -----------------------------------------------------


# construct DiD variables of interest and subset to eligible for treatment
dfdid<- df %>%
  # remove non-partisan/nontechnocrat ministers 
  filter(!(party_min ==0 & techno_min ==0),
         !is.na(party_group)) %>%
  # remove queries submitted to ministries that were technocrat all term
  filter(!(minpartypre2013 == "Non Affiliated")) %>%
  # code treat and post variables
  mutate(post = 1- first_cab, # submitted after cabinet switch
         treat = as.numeric(minpartypost2013 == "Non Affiliated")) # ministry became technocrat


# OLS specifications with and without controls
t2_m1<- (lm(received_resp ~ treat * post, dfdid))
t2_m2<- lm(received_resp ~ treat*post +
             late_in_term + # timing
             sovereignty +large_min +  # ministry
             georeference + casework   + critical + # content
             party_group + asker_leader  + asker_natlist , # asker
           dfdid)

# ses clustered at mincode level
se_t2_m1<- sqrt(diag(vcovCL(t2_m1, cluster = ~ mincode)))
se_t2_m2<- sqrt(diag(vcovCL(t2_m2, cluster = ~ mincode)))


# survival specifications with and without controls
mod<- Surv(dfdid$life_pre, dfdid$received_resp)

t2_m3<- (coxph(mod ~ treat*post, dfdid, cluster = mincode))
t2_m4<- (coxph(mod ~ treat*post +
                 late_in_term + # timing
                 sovereignty +large_min +  # ministry
                 georeference + casework   + critical + # content
                 party_group + asker_leader  + asker_natlist, 
               dfdid,
               cluster = mincode))


## TABLE 2 ##
sink("Table2.tex")
stargazer(t2_m1, t2_m2, 
          t2_m3, t2_m4,
          se = list(se_t2_m1, se_t2_m2, 
                    summary(t2_m3)$coefficients[,4],
                    summary(t2_m4)$coefficients[,4]),
          keep = ":",
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Treat x Post"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(
            c("Controls:", rep( c("", "\\checkmark"),2)),
            c("Model:", rep("OLS", 2), rep("Cox PH", 2))),
          model.names = F)
sink()
####


## full output for appendix

## TABLE A5 ##
sink("TableA5.tex")
stargazer(t2_m1, t2_m2, 
          t2_m3, t2_m4,
          se = list(se_t2_m1, se_t2_m2, 
                    summary(t2_m3)$coefficients[,4],
                    summary(t2_m4)$coefficients[,4]),
          omit = c("mincode*", "party*", "year*", "x"),
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Treat", "Post",
                               "Late in term", 
                               "Ministry of sovereignty","Large ministry",
                               "Georeference", "Casework", "Critical/sensitive",
                               "Asker - leader", "Asker - national list",
                               "Treat x Post"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c( "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(c("Party FE:", rep( c("", "\\checkmark"),2)),
                           c("Model:", rep("OLS", 2), rep("Cox PH", 2))),
          model.names = F,
          no.space  = T)
sink()
####


# Kaplan-Meier plot - MAIN ------------------------------------------------


fit_untreat<- survfit(Surv(life_pre, received_resp) ~ post, filter(dfdid, treat == 0, life_pre >= 0))
untreated<- tidy(fit_untreat) %>% 
  mutate(post = gsub("=", "", str_extract(pattern = "=\\w+", strata)),
         Cabinet = plyr::mapvalues(post, c(0,1), c("2012", "2013")),
         cat = "Untreated") 


fit_treat<- survfit(Surv(life_pre, received_resp) ~ post,  filter(dfdid, treat == 1, life_pre >= 0))
treated<- tidy(fit_treat) %>% 
  mutate(post = gsub("=", "", str_extract(pattern = "=\\w+", strata)),
         Cabinet = plyr::mapvalues(post, c(0,1), c("2012", "2013")),
         cat = "Treated") 

## FIGURE 2 ##
fig2<- bind_rows(untreated, treated) %>%
  mutate(cab = factor(paste0(cat, ", ", Cabinet),
                      levels = c( "Untreated, 2012", "Untreated, 2013","Treated, 2012" ,  "Treated, 2013"  ))) %>%
  ggplot(aes(time, estimate, color = cab, lty = cab)) +
  geom_line() +
  theme_bw() +
  ylim(c(0,1)) +
  xlim(c(0, 1000)) +
  ylab("") +
  scale_color_manual(values = c( "#1F78B4","#1F78B4", "#FC8D62", "#FC8D62"), name = "") +
  scale_linetype_manual(values = c( 1,2,1,2), name = "") +
  xlab("Days after submission") +
  theme(legend.position = "bottom") +
  ylab("Proportion unanswered")

ggsave(fig2, filename = "Figure2.pdf",
       width = 8, height = 5)
####




# DiD analysis - ROBUSTNESS -----------------------------------------------


# more restrictive DiD sample
dfdid_delta<- dfdid %>% 
  # include only ministries that started partisan and changed hands in 2013
  filter(minpartypre2013 != minpartypost2013 & 
           # excludes mins that do not have an analogue in the other period (ie did not exist pre-cabinet shuffle)
           minpartypre2013 != "" &
           minpartypost2013 != "Non Affiliated (PI)")


# OLS models with and without controls
ta7_m1<- (lm(received_resp ~ treat * post, dfdid_delta))

ta7_m2<- lm(received_resp ~ treat*post +
              late_in_term + # timing
              sovereignty +large_min +  # ministry
              georeference + casework   + critical + # content
              party_group + asker_leader  + asker_natlist, # asker
            dfdid_delta)

se_a7_m1<- sqrt(diag(vcovCL(ta7_m1, cluster = ~ mincode)))
se_a7_m2<- sqrt(diag(vcovCL(ta7_m2, cluster = ~ mincode)))


# survival models with and without controls
mod<- Surv(dfdid_delta$life_pre, dfdid_delta$received_resp)

ta7_m3<- (coxph(mod ~ treat*post, dfdid_delta, cluster = mincode))

ta7_m4<- (coxph(mod ~ treat*post +late_in_term + # timing
                  sovereignty +large_min +  # ministry
                  georeference + casework   + critical + # content
                  party_group + asker_leader  + asker_natlist, 
                dfdid_delta, cluster = mincode))

## TABLE A7 ##
sink("TableA7.tex")
stargazer(ta7_m1, ta7_m2, 
          ta7_m3, ta7_m4,
          se = list(se_a7_m1, se_a7_m2, 
                    summary(ta7_m3)$coefficients[,4],
                    summary(ta7_m4)$coefficients[,4]),
          keep = ":",
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Treat x Post"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c("rsq", "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(
            c("Controls:", rep( c("", "\\checkmark"),2)),
            c("Model:", rep("OLS", 2), rep("Cox PH", 2))),
          model.names = F)
sink()
####



# DiD analysis - ASSUMPTIONS ----------------------------------------------


## parallel trends - look at response rate by month for treated/untreated units
a<- dfdid %>% 
  mutate(month_submitted = floor_date(exactpost, unit = "month")) %>%
  group_by(month_submitted, treat, first_cab) %>%
  summarise(rate = mean(received_resp_ever), type = "Any Response",
            sample = "Sample A")
b<- dfdid_delta %>%   
  mutate(month_submitted = floor_date(exactpost, unit = "month")) %>%
  group_by(month_submitted, treat, first_cab) %>%
  summarise(rate = mean(received_resp_ever), type = "Any Response",
            sample = "Sample B")
c<- dfdid %>%
  mutate(month_submitted = floor_date(exactpost, unit = "month")) %>%
  group_by(month_submitted, treat, first_cab) %>%
  summarise(rate = mean(received_resp), type = "Response by Cabinet Enddate",
            sample = "Sample A")
d<- dfdid_delta %>% 
  mutate(month_submitted = floor_date(exactpost, unit = "month")) %>%
  group_by(month_submitted, treat, first_cab) %>%
  summarise(rate = mean(received_resp), type = "Response by Cabinet Enddate",
            sample = "Sample B")

## FIGURE A1 ##
figa1<- bind_rows(a, b, c, d) %>% 
  ggplot(aes(x = month_submitted, y = rate, color = as.factor(treat), fill = as.factor(first_cab))) + 
  geom_point(size = .7) + geom_smooth(span = 1, se = F) + 
  geom_vline(aes(xintercept = minchange))  +
  theme_bw() +
  facet_grid( sample ~ type) +
  xlab("Month Submitted") +
  ylab("Proportion Answered") +  
  scale_color_manual(values = c( "#1F78B4", "#FC8D62"), name = "Treated Ministry",
                     labels = c("No", "Yes")) +
  guides(fill = F) +
  theme(legend.position = "bottom") 

ggsave(figa1, file = "FigureA1.png",
       width = 8, height = 5)
####



## data generating process - prop of queries submitted to treated ministries pre/post switch
dgp1<- dfdid %>%
  mutate(yearmonth = as.numeric(gsub("-", "", substr(exactpost, 1, 7))),
         period = dense_rank(yearmonth),
         qtr = floor_date(exactpost, unit = "month")) %>%
  group_by(qtr) %>%
  summarise(t = mean(treat),
            s = "Sample A")

dgp2<- dfdid_delta %>%
  mutate(yearmonth = as.numeric(gsub("-", "", substr(exactpost, 1, 7))),
         period = dense_rank(yearmonth),
         qtr = floor_date(exactpost, unit = "month"))  %>%
  group_by(qtr) %>%
  summarise(t = mean(treat),
            s = "Sample B")

## FIGURE A2 ##
figa2<- dgp1 %>%
  bind_rows(dgp2) %>%
  ggplot(aes(x = qtr, y = t)) + geom_point() + geom_smooth(span = 1, se = F) + 
  geom_vline(aes(xintercept = minchange)) +
  facet_grid(s ~.) +
  xlab("Month Submitted") +
  ylab("Treated Proportion of Total Queries") +
  theme_bw()

ggsave(figa2, file = "FigureA2.png",
       width = 8, height = 5)
####


## data generating process - query content across treated/untreated,first_caband post switch

# new - critical words
did_content<- dfdid %>%
  mutate(month_submitted = floor_date(exactpost, unit = "month")) %>%
  group_by(month_submitted, treat) %>%
  # monthly averages of casework, region, municipal, and critical/sensitive words
  summarise(Casework = mean(casework),
            Region = mean(reg_ref),
            Municipality = mean(comm_ref),
            critical = mean(critical)) %>%
  mutate(Treat = factor(treat)) %>%
  gather(key, value,Casework:critical) %>%
  mutate(labs = plyr::mapvalues(key, from = c("Casework", "Region", "Municipality", "critical"),
                                to = c("1. Casework", "2. Region", "3. Municipality", "4. Critical/Sensitive"))) 

## FIGURE A3 ##
figa3<- did_content %>%
  ggplot(aes(x = month_submitted, y = value, color = Treat)) + geom_point() + geom_smooth() + 
  facet_wrap(~labs, scales = "free_y") +
  geom_vline(aes(xintercept = minchange)) +
  xlab("Month Submitted") +
  ylab("") +
  scale_color_manual(values = c( "#1F78B4", "#FC8D62")) +
  theme_bw()

ggsave(figa3, file = "FigureA3.png",
       width = 8, height = 5)
####



# observational analysis - ROBUSTNESS -------------------------------------


# OLS models with and without controls + FEs
obs1<- lm(received_resp ~ techno_min, df)

obs2<- lm(received_resp ~ techno_min  +  
            as.factor(year) +late_in_term + # timing
            sovereignty+  large_min +first_cab+  # ministry
            georeference + casework   + critical  + # content
            party_group + asker_leader  + asker_natlist, df)

obs3<- lm(received_resp ~ techno_min  +  mincode +
            as.factor(year) +late_in_term + # timing
            sovereignty + large_min +first_cab+  # ministry
            georeference + casework   + critical  + # content
            party_group + asker_leader  + asker_natlist, df)

# ses at ministry and cabinet level
seobs1<-  sqrt(diag(vcovCL(obs1, cluster = ~ min_period)))
seobs2<-  sqrt(diag(vcovCL(obs2, cluster = ~ min_period)))
seobs3<-  sqrt(diag(vcovCL(obs3, cluster = ~ min_period)))

# survival models with and without controls
obs_cox <- df %>% filter(life_pre >= 0)
mod<- Surv(obs_cox$life_pre, obs_cox$received_resp)

obs4<- (coxph(mod ~ techno_min, 
              cluster = min_period, obs_cox))

obs5<- coxph(mod ~ techno_min +  as.factor(year) +late_in_term + # timing
               sovereignty + large_min +first_cab+  # ministry
               georeference + casework   + critical + # content
               party_group + asker_leader  + asker_natlist , # asker
             cluster = min_period, obs_cox)

obs6<- coxph(mod ~ techno_min + mincode + as.factor(year) +late_in_term + # timing
               sovereignty + large_min +first_cab+  # ministry
               georeference + casework   + critical + # content
               party_group + asker_leader  + asker_natlist, # asker
             cluster = min_period, obs_cox)


## TABLE A8 ##
sink("TableA8.tex")
stargazer(obs1, obs2, obs3, obs4, obs5, obs6,
          se = list(seobs1, seobs2, seobs3,
                    summary(obs4)$coefficients[,4], 
                    summary(obs5)$coefficients[,4],
                    summary(obs6)$coefficients[,4]),
          omit = c("mincode*", "party*", "year*", "x"),
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister",
                               "Late in term", 
                               "Ministry of sovereignty", "Large ministry",
                               "2012 cabinet",
                               "Georeference", "Casework", "Critical",
                               "Asker - leader", "Asker - national list"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c( "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(
            c("Year FE:", rep( c("", "\\checkmark", "\\checkmark"),2)),
            c("Party FE:", rep( c("", "\\checkmark", "\\checkmark"),2)),
            c("Ministry FE:", rep( c("", "", "\\checkmark"),2)),
            c("Model:", rep("OLS", 3), rep("Cox PH", 3))),
          model.names = F,
          no.space = T)
sink()
####


# check whether technocrat effect is the result of co-partisan bias

# OLS with techno and align
al_m1<- lm(received_resp ~ techno_min +  copartisan, df)

al_m2<- lm(received_resp ~ techno_min + copartisan +
             as.factor(year) +late_in_term + # timing
             sovereignty + large_min +first_cab+  # ministry
             georeference + casework   + critical + # content
             soloasker + party_group + asker_leader  + asker_natlist, df)


al_m3<- lm(received_resp ~ techno_min  + copartisan +
             mincode +
             as.factor(year) +late_in_term + # timing
             sovereignty + large_min +first_cab+  # ministry
             georeference + casework   + critical + # content
             soloasker + party_group + asker_leader  + asker_natlist, df)



se_al_m1<-  sqrt(diag(vcovCL(al_m1, cluster = ~ min_period, fix = T)))
se_al_m2<-  sqrt(diag(vcovCL(al_m2, cluster = ~ min_period, fix = T)))
se_al_m3<-  sqrt(diag(vcovCL(al_m3, cluster = ~ min_period, fix = T)))


# survival models
mod<- Surv(obs_cox$life_pre, obs_cox$received_resp)

al_m4<- (coxph(mod ~ techno_min + copartisan,
               cluster = min_period, obs_cox))

al_m5<- coxph(mod ~ techno_min +  copartisan +
                as.factor(year) +late_in_term + # timing
                sovereignty + large_min +first_cab+  # ministry
                georeference + casework   + critical + # content
                soloasker + party_group + asker_leader  + asker_natlist , # asker
              cluster = min_period, obs_cox)

al_m6<- coxph(mod ~ techno_min + copartisan +
                mincode +  as.factor(year) +late_in_term + # timing
                sovereignty + large_min +first_cab+  # ministry
                georeference + casework   + critical + # content
                soloasker + party_group + asker_leader  + asker_natlist, # asker
              cluster = min_period, obs_cox)



## TABLE A10 ##
sink("TableA10.tex")
stargazer(al_m1, al_m2, al_m3, al_m4, al_m5, al_m6,
          se = list(se_al_m1, se_al_m2, se_al_m3,
                    summary(al_m4)$coefficients[,4], 
                    summary(al_m5)$coefficients[,4],
                    summary(al_m6)$coefficients[,4]),
          keep = c("techno*", "copartisan*"),
          dep.var.labels = c("Received Response", "Time to Response (Hazard)"),
          covariate.labels = c("Technocrat Minister", "Co-Partisan Query"),
          
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          float = F,
          omit.stat = c( "f", "ser", "wald", "lr", "logrank", "ll", "max.rsq", "adj.rsq"),
          header = F,
          add.lines = list(
            c("Controls:", rep( c("", "\\checkmark", "\\checkmark"),2)),
            c("Ministry FE:", rep( c("", "", "\\checkmark"),2)),
            c("Model:", rep("OLS", 3), rep("Cox PH", 3))),
          model.names = F,
          no.space = T)
sink()
####


# summary statistics - DATA -----------------------------------------------

# note - the tables in the text were formatted by hand
# the code below outputs the statistics documented in the tables

## TABLE A2 - left panel ##
sink("TableA2_left.tex")
df %>% select(received_resp, techno_min, year, late_in_term,
                  first_cab, sovereignty, large_min, georeference, 
                  casework , critical, 
                  copartisan, asker_leader ,asker_natlist) %>%
  as.data.frame() %>% stargazer()
sink()

## TABLE A2 - right panel ##
sink("TableA2_right.tex")
dfdid %>%
  select(received_resp, techno_min, year, late_in_term,
         first_cab, sovereignty, large_min, georeference, 
         casework , critical, 
         copartisan, asker_leader ,asker_natlist) %>%
  as.data.frame() %>% stargazer()
sink()
####

## TABLE A3 ##
sink("TableA3.tex")
summary(dfregex_main)[[4]][,1:3]
sink()
####

log_close()

